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A systematic asymptotic expansion is developed for the gravitational wave degrees of freedom 
of a class of expanding, vacuum Gowdy cosmological spacetimes. In the wave map description 
of these models, the evolution of the gravitational wave amplitudes defines an orbit in the target 
space. The circumference of this orbit decays to zero indicating that the asymptotic spacetime 
is spatially homogeneous. A prescription is given to identify the asymptotic cosmological model 
for the gravitational wave degrees of freedom using the asymptotic point in the target space. The 
remaining metric function of the asymptotic cosmological model, found by solving the constraints, 
is determined by an effective energy of the gravitational waves rather than from the asymptotic 
point in the target space. 
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I. INTRODUCTION 



Current understanding of cosmology relies heavily on the assumption that the large-scale dynamics of the observable 
universe may be described by an averaged matter distribution evolving in a spatially homogeneous (and isotropic) 
J — ' background. Yet it is clear that the universe is characterized by spatial inhomogeneities that may be large and might 
be dynamically relevant. Some studies have been made to provide a rigorous basis for the validity of the standard 
cosmological approach but the question is not settled 1^, ||, ||, U, |5[ |[ Q . 
^ \ The Gowdy cosmological spacetimes ||] have proven to be valuable theoretical laboratories to explore the behavior 
of spatially inhomogeneous cosmologies. For example, polarized and generic Gowdy models on x i? provided the 
first detailed demonstration of the Belinskii et al (BKL) claim that, in the collapse direction, each spatial point 
— , of a generic cosmological spacetime evolves toward the singularity as a separate universe |l^ . 

' Polarized Gowdy models have also been studied in the expanding direction, both for ||^ and E? spatial 
C\| , topology, where it was argued that the asymptotic behavior could be understood as a stiff fluid of gravitational waves 
■ propagating in a spatially homogeneous background spacetime. Related numerical studies of fluids in plane symmetric 
^ \ expanding cosmologies have also been discussed JTs] , 16 . In this paper, we will argue that this type of asymptotic 



O^- picture, previously obtained for the polarized case, persists for a large class of generic Gowdy models. We shall also 

demonstrate how the background cosmology may be extracted from a numerical simulation of these spacetimes. 
bJO- Ringstrom has recently shown that expanding Gowdy spacetimes fall into two classes with qualitatively different 
behaviors. Only one of these classes is considered here. 



X 

• II. GOWDY MODEL ON x i? 

The vacuum Gowdy cosmology on T'^ x _R is described by the metric [|[ ^ |l3| 

where Q < t < oo, S [0, 27r], {x, y} £ T^, and the gravitational wave amplitudes P, Q and the background variable 
A depend only on 9 and t. If one assumes that P and Q are small, it is easy to see that P is related to the amplitude of 
the + polarization and Q to the x polarization of gravitational waves. (See ||, |l^ for discussion of the interpretation 
of these variables.) Einstein's equations become (with the notation that z,t — dz/dt) 

P,tt + jP,t - P,ee + e^"" {Q,l - Q,?) = (2) 



and 



jQ,t-Q,ee+2P,tQ,t-2P,eQ,e^0 (3) 
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for the wave amplitudes and 



'2 '-e2P( 



= 



(4) 



and 



A,e -2t{P,tP,e 



) = 



(5) 



for the background |25 . Note that Eqs. (Q) and (^) are respectively the Hamiltonian and momentum constraints 
and that A may be constructed after the wave equations (H) and (0) are solved for P and Q. The periodic boundary 
conditions imposed on the metric by spatial topology require M 



A,e d9 = 



(6) 



which may be understood as a restriction on P and Q through Eq. ^ . This condition must be imposed on the initial 
data and is preserved by the evolution. 

Numerical simulation is straightforward since the spiky features which are prominent in collapse |12, |l^ do not 
develop in the expansion direction. However, in contrast to the dominance of local behavior in collapsing models, the 
global behavior is dynamically important so that, at any t, the spatial waveforms must be obtained accurately. 

If we specialize to the case that P, Q, and A depend on t only, Einstein's equations reduce to the spatially 
homogeneous equations (where overdot — d/dt) llj] 



P + 



1 ■ 



and 



„2P/i2 



Q + -Q + 2PQ = 0, 



(7) 
(8) 



X - t[P^ + e^-Pg^) =0 



with the known spatially homogeneous solution ||l7|| 
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(A) 



(11) 



and 



A = Ao + In t 



(12) 



where /x, ^, v, tg, and Aq are constants. For convenience, consider only the case v > as t oo. The metric becomes 
(where some constants have been absorbed and with Q c as t ~* oo) 



ds^ 



(13) 



The contribution due to Q (i.e. the constant c) rotates the coordinate axes with respect to two of the Killing directions. 
In terms of comoving proper time, T, defined by 



(where again some constants have been absorbed), the metric (|l3|) becomes 

ds^ = -dT^ + r2(.^-l)/(.^+3)^^2 ^ j,4(l+.)/(.^+3)(^^ _^ ^j^yY + T^i^-)/(-'+3)dy^. 



(14) 



(15) 
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The exponents of the scale factors T'^' are clearly seen to satisfy the Kasner relations X]fe=i ki = 1 = ^i- We 

shall see later that the presence of spatial inhomogeneities (gravitational waves) will modify this Kasner background 
spacetime. 

As is well-known (see e.g. [0), the Gowdy wave equations (||) and (||) may be regarded as wave maps with the 
target space described by the metric 

dS^ = dP'^ + e^^dQ^. (16) 

A spatially homogeneous solution {P, Q} (t) is represented by a point in the target space while a spatially inhomo- 
geneous solution {P,Q}{t,9) yields a closed curve in the target space. If, asymptotically as t ^ cxd, a spatially 
inhomogeneous solution approaches a spatially homogeneous one, the circumference of the curve will asymptotically 
approach zero and the asymptotic point thus obtained will identify the asymptotic spatially homogeneous solution. 
The symmetries of the target space yield three constants of the motion for the spatially homogeneous equations (1^), 

A = TTqQ - TTp, (17) 
B = TTQ, (18) 

C = TTQ - + 2Qttp, (19) 

where ttp = —tP,t and ttq — —te^^Q,t are the momenta canonically conjugate to P and Q. These constants may be 
related to the constants of the spatially homogeneous solution ( p^ and (|Tl|) through 



Although A, B, and C are not constants of the motion for the wave equations (||) and (|^), their spatial averages A, 
B, and C are constants of the motion. The dynamics of spatially inhomogeneous Gowdy models falls into two classes 
depending on the sign of A^ + BC Jp2t . Here we consider only A^ + BC > 0, the only case allowed for solutions to the 
spatially homogeneous equations (Q) ind (||) and the case assumed in the definition of v in Eq. (|2^). We shall show 
later that knowledge oi A, B, and C allows the asymptotic point in the target space associated with the asymptotic 
spatially homogeneous solution to be identified. 

III. PERTURBATION ANALYSIS 

To understand the asymptotic {t —* oo) behavior of this class of expanding Gowdy models and the dynamical role 
of the gravitational waves assume that P and Q may be expressed as asymptotic expansions 

P^Poit) + epi{e,t) + e^P2{0,t) + (21) 

Q^eqi{e,t) + e^q2i9,t) + ... (22) 

where e is a small parameter and the order term in Q is set equal to zero.. These asymptotic forms are substituted 
into Eqs. (^) and (3|) and terms collected according to their order in powers of e. 

In the development of this perturbation series, it is convenient to insist as we have done here that the order 
solution be polarized — i.e., we require qo{t) ^ ~ qo{t). Then the order equations are just the polarized version of 
Eqs. (^ and (|^) which reduce to 

Po + ^Po = (23) 

with the general solution 



Po = a + vhit 



(24) 
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where a and v are constants. This restriction on the order e solution is not a loss of generality within this class of 
solutions since any solution to Eqs. (0) and (||) with + BC > may be obtained from a polarized solution through 
an SL(2,R) transformation (See, for example, @.) 
To order e^, Eqs. (||) and become 

Pi,tt + 7Pi,t -Pi,ee = >Cppi = 0, (25) 



qi,tt + Q + 7") ~ = -^Q 91 = 0- (26) 

(Note that the second term within the parenthesis in Eq. ( p6| ) is 2po-) The general solutions to these linear equations 
are (with arbitrary constants chosen to eliminate the spatially homogeneous modes of pi and qi) 

N 

Pi^^bn Zo{nt) cos{ne + (j)n), (27) 



N 



(28) 



where n = 1, 2, . . . , A^, a„, &„, 0„, and tpn are constants with Z^{x) any Bessel function of order v |2^. We note that 
the regular and irregular Bessel functions, J^{z) and Y^{z) respectively, have the following large- argument asymptotic 
expansions 



J,(z) = \ — 



TT TT Av — 1 IT TT 

coslz 1> ) smiz V ) 

^ 2 a' 8z ^ 2 4^ 



(29) 



sin(z 



4) + ^ 



■ cos(z - 



1^ 

4' 



(30) 



sm(z- -V- -) 



41^2+3 , TT TT^ 

^^cos(z --.--) 



(31) 



, TT TT, + 
cos(z — —V 7) 

^ 2 4' 8z 



3 . , TT TT 

-sm(z- -V- 



(32) 



where a prime indicates the derivative with respect to the argument. These asymptotic forms for the Bessel functions 
indicate that (1) an arbitrary superposition of Bessel functions of order v and argument nt can be represented 
asymptotically as 



Zy{nt) COs{nt + ifn) 

nnt 



(33) 



where the constant includes any normalization and is an arbitrary phase; (2) as i ^ 00, the leading terms in 
Eqs. ( |27| ) and ( |28| ) of the form Z^{nt) cos{n9 + (5„) have the property that the leading terms of the derivatives with 
respect to t (or, trivially, 6) decay as — i.e., to leading order in inverse powers of t, the derivative acts only on 
the trigonometric functions. 

To order e^, Eqs. (pi) and (0) become 



CpP2 = e^Po (qi,t -qi,g 



(34) 



Cq q2 = 2pi,g qi,e - 2pi,t qi,t ■ 



(35) 
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The order solutions may be expressed as Fourier series (chosen, for convenience, to be the sum of a finite number 
of terms) , 



N 



(36) 



k=l 



N 



92 (, 



,i) =<?f(i) + ^<7r(i)cos(fc0 + Vfe). 



(37) 



fc=i 



The spatially homogeneous modes ^2*^^ q^"^ could, in principle, seriously modify the background spacetime if they are 



(38) 



not small compared to pq (and, in general, qq). Integration over the circle of Eqs. (B4) and (pSD yields 



..(0) 
P2 



■P2 



(0) 



1 

2^ 



if' 



l + 2v 
t 



= \ j {Pi,eqi,e -Pi,tqi,t ) = Sgit), 



(39) 



Wronskian methods [E0| may be used to solve these equations. Note that the Wronskian for Bessel's equation is 



W - M^)Y,:iz) - Y,{z)Jl{z) - — . 



(40) 



If a second order homogeneous ordinary differential equation (ODE) has linearly independent solutions yi{t) and 
2/2(0 ' then the general solution to the inhomogeneous ODE with a source term s{t) is an arbitrary linear combination 
of Di and 2/2 plus the particular solution 



s{t')dt' 



yi{t')v2{t) - V2{t')vi{t) 



(41) 



where the quantity in the denominator is the Wronskian. 

First apply Eq. ( |4l| ) to Eq. (|3^). Linearly independent solutions to the (mathematically) homogeneous equation 
are yi — 1 and y2 = ln< with Wronskian W = 1/t while the source s{t) is just Sp{t). We find 



p(°)=lnt/ Sp{t')t'dt' - I Int' Sp{t')t'dt' 



Partial intgration of the second term cancels the first term on the right hand side (RHS) and leaves 

nt 1 nt' 



pf\t)= j dt'^ J Sp{t")t"dt" 



(42) 



(43) 



Now pi and qi from Eqs. (^^ and (|2^) and their space and time derivatives may be used to construct the source 
terms Sp[t) and Sglt). Note that, from Eq. (p8|). 



N 



qi,t = --qi + t " ^ dnZ±y{nt) cos{n0 + 



(44) 



From Eqs. ( |31| ) and (p2D, it is clear that, as t — > 00, the first term on the RHS may be neglected compared to the 
second. Such terms arising from the time derivative of the factor t^^ in qi will be neglected in the following discussion. 
We find that the right hand sides of Eqs. ( |3^ ) and (35) yield the spacetime dependent source terms 



N N 

Sp{9,t) K. e^" [((i„c?„inmZ^(ni)Z^(mt) cos(n6' + -!/;„) cos(m6' + V'm) 

n— 1 m— 1 

— dndmmnZy(nt)Zc{mt) sin{n0 + ipn) sm{m9 + ipm)] , 



(45) 
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N N 

SgiOjt) « — [bndmnmZQ{nt)Z'^{mt) cos{n9 + 0„) cos{m9 + tpra) 

n—l m— 1 

— bndm'mnZQ{nt)Zy{mt) sm(n9 + tpn) sm{m9 + 4>m)] ■ (46) 

Of course, the products of sines and cosines may be written as cosines of the sums and differences of the arguments. 
Integration over the circle will eliminate the sum of argument terms. We thus obtain the time dependent source terms 



1 ^ 

Sp{t) « 2 5] dln^ e^" [{Z'Ant))' (^.(ni))^] , (47) 



71=1 



1 ^ 

Sgit) « ^ 6„d„n2 cos(V'„ - ^n) [Z^{nt)Zl{nt) - Zoint)Z,{nt)] . (48) 



n=l 



Assume that the mixture of J^, and Y^, is accounted for by an arbitrary phase in the arg ument of the trigonometric 
function in the leading term of the asymptotic expansion. Thus, let Zy{nt) -^2/ imt sin(nt + Cn) and Z'^{nt) k, 
\/2jwni cos{nt + Cn) so that 

N N 

Sp{t) ~^dl^ e^" [cos^ {nt + C„) - siii^ {nt + C„)] ~^dl^ e^" cos(27it + 2C„) . (49) 

n— 1 n=l 

Substitution into Eq. (^) yields 

« E < / ^^'t77 [sin(2nt' + 2C„) - x«] (50) 

where Xn is a constant of integration. 
Asymptotically, 

/■* sin(ai' + 6) , , cos(ai + 5) 

y — f — " — ^ ^''^ 

since the difference between the derivative of the RHS of Eq. dsi] ) and the integrand on the left hand side goes as t^^. 
Using this argument. 



AT 



P2 



(0) 



■ sin(2nt + 2C„) — ^ Ini + ct„ 



47rnt 27r 



(52) 



where cr„ is a constant. It is clear that ^2"^ is indeed a perturbation of po = vlnt since an arbitrary amount of the 
homogeneous solution of Eq. (|38|) may be added to eliminate the logarithmic and constant terms in Eq. (m) . 



Now solve Eq. (|39|) for ^2°'' • Here the linearly independent solutions arc yi = 1 and y2 = t ^" with the Wronskian 
W — —2vt^'^'"^^. Substitution into Eq. ( p] ) shows that the particular solution to Eq. ( |39| ) is 

(i) = -^t-'^ f SQ{ty'"+'dt' + SQ{t')t'dt' . (53) 

Taking the time derivative and then integrating shows that Eq. ( p3[ ) is equivalent to 

gf(t)= f dt'{t'y^^-^ f dt"{t" f^+^SQ{t"). (54) 



Just as before, we evaluate the asymptotic form of SQ{t). With the previous asymptotic form for Z„ and its derivative 
and Zo{nt) « ^Ij-Knt sin(ni + ^„) and Z'^ipi) k, 2 /Tint cos{nt + ^„), 

N 

^Qi^) ~ E [cos(n-t + g„) cos(»i + Cn) - sm{nt + g„) sin(n^ + Cn)] 

n—l 

W 

w E ^ndn ^^^^^ cos(2nt + C« + Cn) ■ (55) 
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After partial integration, the integral over t" in (^J) yields as the dominant term 

J2 Kdn j dt' {t')-"'' sin(2nt' + C„ + Cn) ■ (56) 



N 

(0) 
92 



A multiple of the homogeneous solution has been added to eliminate constants of integration. The same argument 
[Eq. (pl|)] concerning asymptotics of the integrals is valid here. Essentially, the leading term in Eq. (^^ is obtained 
by integrating only the trigonomentric term. Thus we obtain 



« ^6„d„^— ^cos(2nt + C„+^„). (57) 



Again, this quantity is seen to be a decaying perturbation on go's constant value. 

The solutions ^2"^ and q^'^ decay as and t~'"~^ respectively — i.e. faster than pi and qi — and are perturbations 
of po and qg. [For a general (unpolarized) order e" solution, as t 00, pq vlnt and qo —> c, a, constant.] The next 

relevant order for the spatially homogeneous mode, p^^'^ and will have source terms containing (e.g.) e^" {q')iPi 
that will decay as t~^/^, with every other aspect of the argument the same. 

Now consider the spatially dependent mode amplitudes pj'^'' and gj'^^ for fc > 1. Integrate cos{k6 + ipk) Sp{0, t) over 
the circle to yield the fcth component 

N 

Sp\t) « ^[dndik-nin-lk - n\Z'^{nt)Z'^{\k - n\t) 

n=l 

-dnd\k-n\n\k - n\Zy{nt)Zy{\k - n\t)] 

N 

+ [dndk+nn{k + n)Z[,{nt)Z'^{{k + n)t) 

-dndk+rMk + n)Zy{nt)Zy{{k + n)t)] . (58) 
We again construct the particular solution pj''^ using ( ^ ) where yi — Jo{kt), y2 — Fo(fci), and W — 2/{Trkt) . Then 

pi''^ = Yoikt) f sl^\t')Mkt'f-^dt' ~ Mkt) f S^'^\t')Y^{kt')^dt' . (59) 



We can regard the asymptotic expansions of and to be those given in ( 29 ) and (^ while Zi, and its derivative 
have an arbitrary phase. 
Asymptotically, 

Sp\t) » dnd\k-n\ \k - n\--cos[{n + \k - n\)t + 2^^] \ 

n=l 



^ r 1 

\ dndk+n{k + n)-c 



cos[{2n + k)t + 2xpn\ \ ■ (60) 



Rather than complete the evaluation of Eq. ( pOj ) , notice that power counting gives the dominant asymptotic behavior 
as . The same type of asymptotic argument about the integral of a trigonometric function times a power law that 

(k) 

was used previously implies that P2 should be an oscillatory function divided by t as expected. The argument for 
q^^ should be identical except that appropriate powers of t^"" must be carried along. 

We note that in general, the sums and products of trigonometric functions will yield trigonometric functions. (Any 
potentially troublesome constant term has already been considered in the sources for p^2^ and q'2^ ■) As long as 

there are trigonometric functions, the leading order as t ^ 00 of time derivatives of p^2^ and q"2^ will arise from 
the derivative acting only on the trigonometric term and not on the power law. As we have previously argued, the 

ik') (k^ 

corresponding statement is also true for integrals over time. Thus we expect ^^^id q^ to have the forms 

^ (M) (.) ^ Q^\t.e) 

P2 ^ ■ I2 lOlj 
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where and are oscillatory functions of 9 and t. 

Continuing these arguments yields that at each order we obtain for p„ and g„ functions which decay faster than 
those of the previous order by ^/t as t ^ oo. This behavior arises because, in general, p„ and g„ for n> 2 will be 
the solution to an inhomogeneous wave equation with source terms of the form ("g" j)*^ i"p" lY where kj + Ir = n and 
" " indicates that space and time derivatives may act on qj and pi. (Asymptotically, these derivatives act only on 
the trigonometric functions so that power counting yields the correct power law time dependence. Note that q„ will 
always have an extra factor of in its power law time dependence.) The fundamental quantities in the source terms 
are pi and qi (e.g. appearing quadratically in the source terms for p2 and 52) which decay as i"-'^/^. In addition, from 
Eq. (|6l|), all orders (fc) in the Fourier series for p„ and qn for n > 1) will have the same power law dependence on t. 
Substitution of Eq. ( |6l| ) in the Fourier series (|3^) , ( |37| ) will thus yield 

^ V 

P^a + vlnt + Y,^, (62) 
i=i 

Q-^ + Et^ (63) 



where 

,(fc) 



= J2rpit) cos{ke + ^k), (64) 



fc 



In the rest of this paper, we shall consider evidence for the eventual dominance of the asymptotic behavior given by 
Eqs. (H) and M). 



IV. DECAY TO A SPATIALLY HOMOGENEOUS COSMOLOGY 



The asymptotic forms ( |62| ) and (|6^) suggest that expanding Gowdy spacetimes approach the spatially homogeneous 
cosmology described by (po, go), the solutions ( |lO| ) and ( pT| ) (in the general case). Recall that these solutions are 
completely determined by the constants v, fj,, ^, and to and that the first three of these are related to the constants 
A, B, and C of the spatially homogeneous solution defined in Eqs. (|l3)-(|9l) through Eq. (|^). While A, B, and C 
are not constants of the motion for the wave equations (||), (||), their integrals over the circle, 

A^^<l)Ade, B=^j>Bde, C=^<j)Cd0, (66) 

are constants of the motion. These are the well-known conserved quantities associated to the Killing symmetries of 
the target space of a harmonic map — in this case the hyperbolic plane. If P and Q evolve to po and qo respectively as 
t — * 00, then they will approach the one parameter family of solutions such that A, B, C fix v, p, ^ through Eq. ( [20| ) 
with A replaced by A, etc. The constants A, B, and C determine a point in the target space. For any solution to the 
spatially homogeneous equations (Q) and (|^), + BC > 0. Since the perturbation analysis considered here expands 
about a solution to Eqs. (0) and (g), the signature + BC > is built in. However, the opposite sign is allowed 
for spatially inhomogeneous Gowdy models and has been shown by Ringstrom to yield a qualitatively different 
asymptotic state. 

A, B, and C may be constructed perturbatively using and (H). If we consider such an expansion for any of the 
constants, say B = bnS", we can show by interchanging derivatives with respect to t and e that each coefficient 6„ 
in B must be separately constant. However, in general, 6„ will have a nontrivial power law dependence on t for n > 0. 
For bn for n > to be constant, the coefficient of the power law must vanish to yield B = bo, the value obtained from 
the homogeneous solution po, qg. The same argument may be used for A and C. Furthermore, in general, at least for 
the first few terms in the asymptotic expansion, the formal time dependence may be a growing power law. We again 
consider B as an explicit example. 

N 

B = Y,bne'' ^ -— j dOte'P- [go e" + qi,t e + (2pigi,t +92,* )e' + • ■ ■] • (67) 

n— 
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In the perturbation analysis of the previous section, we have assumed qq = = qq. Thus the order e° term vanishes. 
This yields the expected value for polarized solutions, B = 0. The order e term vanishes because the spatial average 
of qi is zero. Note, however, that the integrand has the overall asymptotic time dependence of ^"+1/2 and thus will 
grow for V > —1/2. The order term, &2, in ( |67|) is 

b2 = ^ j^d0t'+^"i2piqi,t +92,*) 

= t'+^^(4"^ + ^ (fd9 2p,q,,t). (68) 



2tt . 

Unless V < —1/2, the overall t dependence indicates a growing term. However, from (|5^) 



.(0) ^ (^)-2.-i f dt'HT'-^SQit'). (69) 



Substition of the source term written as the integral over the circle of the right hand side of Eq. (35) yields 



'?2°^ = 7TT2^ fdt'^j>de{t'Y+^-{p,,equ9-Pi,t'qi,t'). (70) 
Partially integrating the first term on the RHS with respect to 9 and second with respect to t' yields 



=-2^ /'^^fl?!'* / dt'^jdOp, \^{t'f+^^q,,,,-±[[tT^\,,,,] \ . (71) 

The expression inside the braces is just the wave equation ( |26| ) for qi and thus vanishes to give 

+ j de2piqi,t^Q. (72) 

Note that the values of A, B, and C are determined from the order solution. The order e" terms for n > 1 vanish 
due to the wave equations, even though the integrands A, B, and C may grow in t. 

To identify completely the background spacetime also requires the construction of A through Eq. (^ . We now show 
that the presence of gravitational waves changes the time dependence of the spatially homogeneous mode of A from 
that of the spatially homogeneous solution (|l2|) (see [^). Recall that if, asymptotically, 

P^vlnt + (73) 



then the asymptotic forms of their time derivatives are dominated as i — > oo by 

rather than by po — v/t where Vi^t and Qi,t are oscillatory in time. Similarly, the dominant terms in the spatial 
derivatives are 

P,e«^ , Q,e«-j^. (76) 

Eq. (^ then becomes 

X,t~Vlt +Qlt +Vle +Qle. (77) 

To compute A, the spatial average of A, the bounded, oscillatory functions on the RHS of Eq. (|77| ) may be replaced 
by a constant average value, ^, plus fluctuations to yield 

A « a (78) 
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in contrast to the logarithmic dependence on t in Eq. (|12D. Note that ( depends on the details of the first order 
fluctuations rather than on the constants of the spatially homogeneous solution. 
Given Eq. ([78|), the metric becomes, as t oo, 



ds^ = e'-'lH-^l^ (- de + dd^) + {dx + cdyf + t^-^dy". (79) 
We again must transform to comoving proper time defined by 

dr = e^/2-i/4int^i_ (80) 



An approximate form of the metric ( |79| ) may be found. As < — )■ oo, expressions of interest may be expanded in the 
small quantity (lnt)/t. Thus Eq. (|80| ) becomes 



Keeping only the dominant term yields 



or, neglecting some constants. 



dr«e'^*/* (^1- ^Int^ dt. (81) 



T « ^e^*/4 (82) 



t«^lnT. (83) 



This, in turn, gives the approximate metric (again absorbing constants as needed into redefinitions of the spatial 
variables) 

da^ « -dT^ + T'^d¥ + (lnr)i+"(d£ + cdyf + (lnT)i-"dy2 (84) 

where the tilde indicates rescaling and we have used k, (16/C^) e''*/'^ t^^l'^ . Now consider this approximate metric 
at some time % such that lni(T,) << i(Tj). Then, for T >T^,T ^T^ + (T - T,) = T, + AT where we shall assume 



that AT << Ti. Then InT « InT; + AT/T^. As long as AT/Ti « InT^, the approximate metric (|84|) becomes 

ds^ K -dT^ + T^dP + {liiT,f+''{di + ddyf + {\nTif-'"df. (85) 

If the spatial variables (and c) are defined once again to absorb powers of InT^, the approximate metric ( p5| ) may 
be identified as the cosmological sector of flat spacetime in Rindler coordinates. A more detailed discussion of the 
identification of the limiting metric is given for the polarized case in |9j. 

In the following section, we shall demonstrate numerically the decay of the Gowdy solution to this background 
spacetime. First we shall consider two quantities which require the presence of the spatially inhomogeneous modes of 
P and Q. Both quantities may be understood by starting from the Hamiltonian density 

W = ^ (4 + e-'^-l) + \ {P,l + e''Q.l ) (86) 
whose variation yields Eqs. (|^) and (||). The quantity 

C = fd0{P,l+e'^Q,lf' (87) 

is the orbit circumference in the target space. If C — *■ as i — > oo, the solution approaches a spatially homogeneous 
one. Substitution of the asymptotic forms (|6^), (|6^) into (|8^) shows that we expect 



,83, 

where p(t) is a bounded oscillatory function which is, asymptotically, the spatially homogeneous mode of 
(7^1,9 + Qi,l y^^. Another function expected to decay to zero if the asymptotic spacetime is spatially homogeneous 
is the generalized force 

T deiyn-vny^ (89) 
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where V = {S/5P, e ^6/5Q) is the gradient in the target space. We find 



de 



1/2 



(90) 



which has the asymptotic form 



(91) 



with ^/;(t) a bounded oscillatory function. 

Assuming that C and T decay, the resultant asymptotic point in the target space is one of the spatially homogeneous 
solutions characterized by a particular set of w, /i, ^ obtained from A, C through (|20|). The remaining parameter 
to represents a rescaling of the time t. It may be determined by fitting data obtained from numerical simulations of 
the full Einstein equations for the model. 



V. NUMERICAL SIMULATIONS 



We use numerical simulation to provide evidence in support of our previous analysis of the asymptotic behavior of 
the Gowdy spacetimes. As pointed out previously, the lack of spiky features in the expanding spacetimes means that 
almost any numerical technique can evolve the spacetime. Here we use an iterative Cranck-Nicholson algorithm [^ . 
Many of our conclusions will be based on the spatially homogeneous mode (i.e. the average over the circle) of various 
quantities. Since most of these quantities are oscillatory in both space and time with large amplitude, accurate spatial 
waveforms are needed to generate the correct average. 

The quality of the numerics in the simulation may be found by computing the conserved quantities j4, i3, and C 
These are shown for representative initial data in Fig. |l|. Note that A, B, and C have been singled out from among 
all combinations of constants of the motion of the spatially homogeneous solution because their time derivatives are 
also total derivatives with respect to 9. While v, fj,, and ^ are constants of the spatially homogeneous solution, their 
integrals over the circle in the general case are not constants in time. Fig. |2| shows one of the conserved quantities 
(viz. A) on a much finer scale so that deviations from constant value are visible. These deviations display second order 
convergence to zero (at least during most of the simulation). In contrast. Fig. || shows the integrands of A, B, and 
C at a representative spatial point. These oscillate wildly in time. The observed constancy of A, B, and C indicates 
that the simulation is reliable. 

For a reliable simulation it is then possible to compute C and J- using Eqs. ( ^7| ) and ( po|) respectively. These are 
shown in Fig. ^ along with power law fits to the data. A closeup of the oscillatory functions p{t) and ■0(t) (see Fig. ^b) 
argues that the larger discrepancy in between the expected t~^/^ and the observed decay is due to the greater 
difficulty in determining the centroid of the pattern. Similar decay as t^^^^ is shown for other representative initial 
data in Fig. |[ 

If these solutions indeed decay to a spatially homogeneous background spacetime, it now becomes possible to 
identify the precise asymptotic behavior of P and Q. Since P and Q are periodic functions on [0, 27r], it is always 
possible to compute the spatially homogeneous mode of the appropriate Fourier series. This mode for P and Q (say 
p and q) contains the e° solution as well as the spatially homogeneous modes of the order e" solutions for n > 2. 
However, if the asymptotic forms given in Eqs. ( |6^ ) and ( |63| ) indeed describe the solution, we require (as argued in 
Section III) that the higher order contributions to p and q should be oscillatory (and decaying) and should oscillate 
about the spatially homogeneous solutions pa and qo respectively. Evidence for such behavior is the following: The 
constants A, _B, C may be computed from the initial data and used to construct /i, ^, and v. This yields po and go 
as functions of t/to where to is as yet unknown. Fig. ^ shows that if to is chosen to yield a best fit to (say) p late in 
the simulation, it also fits the correct average time dependence of p and yields the correct asymptotic behavior for q. 
The fact that we can fit two functions with one adjustable parameter provides strong support that the computed /z, 

V, and the fit to define the correct asymptotic P and Q. 

Note that the time averaged values of p and q depart from po and qo early in the simulation. This means that, 
although /i, V, and ^ are fixed by the initial data, to does not have the value one would calculate from the initial data. 
This means that nonlinear interactions which arc important early in the simulation influence (i.e. renormalizc) the 
value of to and thus the background spacetime. 

Also note that it is possible to evaluate to without fitting by constructing p and q defined as the solutions ( [lO| ) and 
( |ll| ) where fi, and v are those found from the initial data while to is evaluated from p, q at that time step. This 
procedure yields io(t) which converges to io as t — > cx). Fig. ^ shows p, g, and to compared to po, qo, and to- 
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Fig. H shows the spatial average of A vs. t for a representative simulation. The time dependence is clearly linear 
(and definitely not logarithmic). Also shown is the spatial average of the energy-like terms, £ in Eq. (J^, which are the 
"source" for A. Detailed examination shows that the fluctuations in the source terms decay (as 1/t) to the appropriate 
constant value. 

The analysis of Sections III and ^ and the numerical simulations have provided support for an asymptotic state 
described by gravitational waves of decaying amplitude propagating in a spatially homogeneous cosmology. Figures 
p[-p^ illustrate the mechanism yielding this asymptotic state. In particular, the source terms Sp(t) and SQ{t) must 
be driven to small amplitude if they are not small initially. Note that the exponential dependence on P in Sp could 
be catastrophic if, for example, P were to depend linearly on t. In fact, such linear behavior of P is seen very early in 
the simulation in Fig. |ll|. There, the nonlinear terms then act as regulators (as described for the collapsing direction 
in p^) to suppress the linear behavior. 

The consistency between the asymptotic behaviors of P, Q, and A found in the numerical simulations and those 
expected from the analysis in Sections III and 



terms Sp {t) and Sq (t) on po {t) and qo (?) must 



[V suggests that the behavior found there for the nonlinear source 



DC correct — i.e. these terms decay as 1/t or faster yielding decaying 



perturbations pn"^ and qii^ for n > 2 of po and qo respectively. It is interesting to examine direct evaluation of the 
relevant terms during the simulation. Figures |9| and |0| illustrate respectively the component terms of Sp and Sq 
(along with the corresponding P and Q) for the entire simulation starting from the initial data of Fig. nl Early in 
the simulation, nonlinear effects dominate. The term in Sp, Vi — § d9e^^ Q,^, in fact behaves as e~ (i.e. it is 
dominated by the constant of motion B) as is seen to be the case in Fig. ^ and, for more complicated initial data, 
in Fig. |ll|. Figures Eo| an d |lp illustrate the complicated evolution of the two terms in Sq which, in the end, are seen 
to decay. Figures |13| and |14| show the evolution of Sp and Sq early and late in the simulation (for the complicated 
initial data). Nonlinear effects are clearly significant early in the simulation. Late in the simulation, however, Sp and 
Sq have decayed to small amplitudes with even smaller (by several orders of magnitude) time averages (see Fig. p^ . 
This oscillatory (in time) behavior of the source terms means that any (positive or negative) growth in P tends to 
produce a correction to reduce the magnitude of the growing term. Thus the time averages of these source terms 
represent their actual, negligible influence on the dynamics. 



VI. CONCLUSIONS 



We have shown here that the asymptotic state of a class of expanding Gowdy models may be described quantitatively 
as decaying amplitude gravitational waves propagating in a spatially homogeneous background spacetime. Quantities 
such as the orbit circumference in the target space, C, and generalized force, T , are seen to decay to the background 
value of zero as as predicted by perturbative methods. (After these results were obtained, Ringstrom was able 

to prove the form (^) for C for "small" initial data |^.) In addition, we have provided a prescription for construction 
of the asymptotic background spacetime. 

Open questions concern the generality of the results discussed here. While the asymptotic behavior (^2|), (|6^) has 
been found for many sets of initial data for this class of Gowdy spacetimes, there is no guarantee that the parameter 
space of such data has been fully explored. (In this context, Ringstrom's result for the decay of the gravitational wave 
amplitudes is quite encouraging.) One hopes that this will be settled with mathematical tools. 

A more significant issue is whether or not other classes of spatially inhomogeneous cosmological spacetimes may 
be described in terms of matter and gravitational radiation with decreasing influence on a background spatially 
homogeneous cosmology. Most cosmological studies assume that the present universe may be described this way. 
Within the Gowdy models considered here, this type of description arises after a nonlinear regime which has an 
influence on the final state of the system. 

Investigations of more general models are in progress. We note that, in contrast to the collapsing case where most 
types of "matter" are dynamically unimportant, every modification of the nature of the expanding model and its 
matter content affects the dynamics. Finally, we remark that a modification of the techniques used here have been 
applied to vacuum, general symmetric spacetimes (see for a description of these spacetimes). Results will be 



discussed elsewhere |24 
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FIG. 1: Constants of the motion A, B, and C in a representative expanding Gowdy spacetime. Unless otherwise noted, all 
figures are constructed from a simulation with P = 1, P,t = 2 + 2cos^, Q = cosO, Q,t = 0.3, and A = 0.6 cos^ at t = 1. 
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FIG. 2: Fluctuations of A. On a fine numerical scale, numerical fiuctuations in A (which should be strictly constant) are 
visible. In the main graph, the three curves show (from top to bottom) fluctuations with spatial resolutions of 4096, 2048, and 
1024 points respectively. Offsets have been added to the lower resolution values to allow display on the same graph. The inset 
shows the mean values where the spatial resolution increases from top to bottom over the same range in time. It is clear that 
numerical errors in the constant A are converging to zero and that the mean value is also converging. 
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FIG. 3: Integrands for A, B, and C at a representative spatial point. While the constants A, B, and C are preserved in the 
numerical simulation, the integrands A, B, and C fluctuate wildly. Here they are rescaled to allow display on a single graph. 
The horizontal line represents the average values. 




FIG. 4: Evidence for approach to a spatially homogeneous background cosmology. The circumference C and force are 
shown. The dashed and solid thick lines axe power-law fits with the indicated exponents. The region inside the rectangle is 
shown in the inset. The lower curve (solid line or line with circles) represents C and the upper curve (dashed line or line with 
squares) represents J^. 




FIG. 5: Decay of C (solid lines) and (dashed lines) for three different expanding Gowdy models. The respective power-law 
fits are shown at the top of the graph. 
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FIG. 6: Predicting the asymptotic spatially homogeneous model, (a) The values of the spatially homogeneous modes, p of 
P (solid line with squares) and q of Q (dashed line with circles), are compared to po and qo calculated from /i, ^, v, and the 
best fit to to ■ (b) Rescaling the time associated with given values of po and go until po matches the center of p. This occurs if 
to = 5.837. (Note that the value of to will change slightly as the match is made for increasing t, eventually reaching the "true" 
value at t = oo.) (c) A comparison between q (dashed) and qo (solid) late in the simulation. Note that q is decaying to qo and 
that the difference is only a few parts in 10*. 
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FIG. 8: Linear growth of A. The dots represent every 20th point in the spatial average of A as computed using Eq. in a 
numerical simulation. The solid line is A = 23.00 + 4.4969 1. Also shown is the energy-like spatial average of the source terms 
for A,t, which, while fluctuating, asymptote to a constant. Note the small range on the scale for E. The thick solid line is 
E = 4.4969 + 5.6388^. 
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FIG. 9: Nonlinear fluctuations and the spatial average of P. The nonlinear terms Vi = ^d9e^^Q,l (thin solid line) and 
V2 = §d9e^^Q,l (dashed line) are shown on a log- log plot (using the vertical axis labeled NL). The logarithmic time axis is 
used to emphasize the behavior during the early stages of the simulation when nonlinear effects arc significant. For comparison, 
the spatial average of P is also shown (thick solid line). Note that Vi increases and V2 decreases as P decreases as expected if 
these terms are acting to regulate P. 




FIG. 10: Nonlinear fluctuations and the spatial average of Q. The nonlinear terms Wi = 2\fddP,tQ,t \ (thin solid line) 
and W2 = 2| ^ d9 P,g Q,g \ (dashed line) are shown on a log-log plot (using the vertical axis labeled NL). For comparison, the 
spatial average of Q is also shown (thick solid line). 
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FIG. 11: Nonlinear terms and the spatial average ol P early in the simulation. The nonlinear terms V\ (thin solid line) and V2 
(dashed line) (see Fig. ^ are shown for 1 < i < 6 (using the vertical axis labeled NL) . For comparison, the spatial average of P is 
also shown (thick solid line). For this and subsequent figures, the initial data are P = 0, P,t = 2 + 2 cos 6 + 0.5 cos(26) + 4 cos(50), 
Q — cos9 + 2cos(36') + cos(7&), Q,t = 0.3. Note that, as P becomes increasingly negative early in the simulation, Vi increases 
while V2 decreases supporting the contention that the former evolves as e~^^ . Also note the qualitative similarities to Fig. ^ 
despite quite different initial data. 




FIG. 12; Nonlinear terms and the spatial average of Q early in the simulation. The nonlinear terms Wi (thin solid line) and 
W2 (dashed line) (see Fig. |l^) are shown for 1 < t < 6 (using the vertical axis labeled NL). For comparison, the spatial average 
of Q is also shown (thick solid line). The horizontal line indicates zero for the nonlinear terms. Note that the difference in 
initial data leads to qualitative differences in evolution in Figs, nd and [13. 
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FIG. 13: Nonlinear source terms Sp(t) (solid line) and SQ{t) (dashed line) early in the simulation. Note how quickly these 
terms approach a state where they oscillate more or less symmetrically about zero. 
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FIG. 14: Nonlinear source terms Sp{t) — Vi — V2 (solid line) and SQ{t) = W2 — Wi (dashed lino) late in the simulation. Note 
the small amplitude of these terms and that their time averages are very close to zero (9 x 10"* and 2 x 10~^ respectively). 



